function [beta,R2] = reg_Y_Nz_on_X_Nz(Y_Nz,X_Nz,prob_Nz)

    mean_X = sum(prob_Nz(:).*X_Nz(:));
    var_X = sum(prob_Nz(:).*(X_Nz(:).^2)) - mean_X^2;
    mean_Y = sum(prob_Nz(:).*Y_Nz(:));
    var_Y = sum(prob_Nz(:).*(Y_Nz(:).^2)) - mean_Y^2;
    mean_XY = sum(prob_Nz(:).*X_Nz(:).*Y_Nz(:));

    beta = (mean_XY - mean_X*mean_Y)/var_X;

    R2 = beta^2*var_X/var_Y;

end